################################################################################
############################    Plot for reliability ###########################
################################################################################
dat <- read.csv("PGY5_4.csv")%>%
  select(-c("rater","PGY","item","X"))
formula.dat <- "score ~ (1 | subjectID) + (1 | procID)"
g1 <- gstudy(data = dat, formula = formula.dat, colname.strata = "subset",
             colname.objects = "subjectID")

############################ PGY2 #########################################
n <- 1:20
#relative error variance
PGY2_c <- read.csv("PGY2_c.csv")
rel_err_var <- PGY2_c[9,4]/n
#absolute error variance
abs_err_var <- PGY2_c[9,3]/n +PGY2_c[9,4]/n
#calculate generalization coefficient
gen_coef <- PGY2_c[9,2]/(PGY2_c[9,2]+ rel_err_var)
#calculate dependability coefficient
dep_coef <- PGY2_c[9,2]/(PGY2_c[9,2]+ abs_err_var)
#calculate SEM
SEM <- sqrt(abs_err_var)

label <- "phi"
data <- data.frame(n=n,dependability=dep_coef,SEM=SEM) %>%
  mutate(label=label)


p_PGY2 <- ggplot(data,aes(x=n,y=dependability)) +
  geom_text(aes(n,dependability,label=label), parse=TRUE) +
  geom_hline(yintercept=0.8, linetype="dashed", color = "red") +
  scale_y_continuous(breaks = c(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0)) +
  xlab("Number of Procedures") + 
  ylab("Reliability Coefficient") + 
  labs(title = "PGY 2", tag = "a") +
  theme_bw()+
  theme(axis.line = element_line(colour = "black"),
        panel.background = element_blank(),
        panel.border = element_blank(),
        legend.title = element_blank())+
  theme(plot.title = element_text(hjust = 0.5,size=10),
        plot.tag = element_text())

############################ PGY3 #########################################
n <- 1:20
#relative error variance
PGY3_c <- read.csv("PGY3_c.csv")
rel_err_var <- PGY3_c[19,4]/n
#absolute error variance
abs_err_var <- PGY3_c[19,3]/n +PGY3_c[19,4]/n
#calculate generalization coefficient
gen_coef <- PGY3_c[19,2]/(PGY3_c[19,2]+ rel_err_var)
#calculate dependability coefficient
dep_coef <- PGY3_c[19,2]/(PGY3_c[19,2]+ abs_err_var)
#label <- "rho"
label <- "phi"
data <- data.frame(n=n,dependability=dep_coef) %>%
  mutate(label=label)


p_PGY3 <- ggplot(data, aes(n=n,y=generalizability)) +
  geom_text(aes(n,dependability,label=label), parse=TRUE) +
  geom_hline(yintercept=0.8, linetype="dashed", color = "red") +
  scale_y_continuous(breaks = c(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0)) +
  xlab("Number of Procedures") + 
  ylab("Reliability Coefficient") + 
  labs(title = "PGY 3", tag = "b") +
  theme_bw()+
  theme(axis.line = element_line(colour = "black"),
        panel.background = element_blank(),
        panel.border = element_blank(),
        legend.title = element_blank())+
  theme(plot.title = element_text(hjust = 0.5,size=10),
        plot.tag = element_text())

############################ PGY4 #########################################
n <- 1:20
#relative error variance
PGY4_c <- read.csv("PGY4_c.csv")
rel_err_var <- PGY4_c[19,4]/n
#absolute error variance
abs_err_var <- PGY4_c[19,3]/n +PGY4_c[19,4]/n
#calculate generalization coefficient
gen_coef <-
  PGY4_c[19,2]/(PGY4_c[19,2]+ rel_err_var)
#calculate dependability coefficient
dep_coef <-
  PGY4_c[19,2]/(PGY4_c[19,2]+ abs_err_var)
#label1 <- "rho"
label1 <- "phi"
data <- data.frame(n=n,dependability=dep_coef) %>%
  mutate(label=label1)


p_PGY4 <- ggplot(data, aes(n=n,y=dependability)) +
  geom_text(aes(n,dependability,label=label), parse=TRUE) +
  geom_hline(yintercept=0.8, linetype="dashed", color = "red") +
  scale_y_continuous(breaks = c(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0)) +
  xlab("Number of Procedures") + 
  ylab("Reliability Coefficient") + 
  labs(title = "PGY 4", tag = "c") +
  theme_bw()+
  theme(axis.line = element_line(colour = "black"),
        panel.background = element_blank(),
        panel.border = element_blank(),
        legend.title = element_blank())+
  theme(plot.title = element_text(hjust = 0.5,size=10),
        plot.tag = element_text())

############################ PGY5 #########################################
n <- 1:20
#relative error variance
PGY5_c <- read.csv("PGY5_c.csv")
rel_err_var <- PGY5_c[33,4]/n
#absolute error variance
abs_err_var <- PGY5_c[33,3]/n +PGY5_c[33,4]/n
#calculate generalization coefficient
gen_coef <-
  PGY5_c[33,2]/(PGY5_c[33,2]+ rel_err_var)
#calculate dependability coefficient
dep_coef <-
  PGY5_c[33,2]/(PGY5_c[33,2]+ abs_err_var)
#label1 <- "rho"
label1 <- "phi"
data <- data.frame(n=n,dependability=dep_coef) %>%
  mutate(label=label1)
#The symbol 𝜌 represents the generalization

p_PGY5 <- ggplot(data, aes(n=n,y=dependability)) +
  geom_text(aes(n,dependability,label=label), parse=TRUE) +
  geom_hline(yintercept=0.8, linetype="dashed", color = "red") +
  scale_y_continuous(breaks = c(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0)) +
  xlab("Number of Procedures") + 
  ylab("Reliability Coefficient") + 
  labs(title = "PGY 5", tag = "d") +
  theme_bw()+
  theme(axis.line = element_line(colour = "black"),
        panel.background = element_blank(),
        panel.border = element_blank(),
        legend.title = element_blank())+
  theme(plot.title = element_text(hjust = 0.5,size=10),
        plot.tag = element_text())

tiff("Generalizability.tiff", unit="in",width=7, height=5, res=500)
ggarrange(p_PGY2,p_PGY3,p_PGY4,p_PGY5)
dev.off()




################################################################################
############################    Plot for SEM  ##################################
################################################################################

############################ PGY2 #########################################
n <- 1:20
#relative error variance
PGY2_c <- read.csv("PGY2_c.csv")
rel_err_var <- PGY2_c[9,4]/n
#absolute error variance
abs_err_var <- PGY2_c[9,3]/n +PGY2_c[9,4]/n
#calculate generalization coefficient
gen_coef <- PGY2_c[9,2]/(PGY2_c[9,2]+ rel_err_var)
#calculate dependability coefficient
dep_coef <- PGY2_c[9,2]/(PGY2_c[9,2]+ abs_err_var)
#calculate SEM
SEM <- sqrt(abs_err_var)

label <- "Delta"
data <- data.frame(n=n,SEM=SEM) %>%
  mutate(label=label)


p_PGY2 <- ggplot(data,aes(x=n,y=SEM)) +
  geom_text(aes(n,SEM,label=label), parse=TRUE) +
  scale_y_continuous(breaks = c(0,0.2,0.4,0.6,0.8,1,1.2,1.4,1.6))+
  xlab("Number of Procedures") + 
  ylab("Standard Error of Measurement") + 
  labs(title = "PGY 2", tag = "a") +
  theme_bw()+
  theme(axis.line = element_line(colour = "black"),
        panel.background = element_blank(),
        panel.border = element_blank(),
        legend.title = element_blank())+
  theme(plot.title = element_text(hjust = 0.5,size=10),
        plot.tag = element_text())

############################ PGY3 #########################################
n <- 1:20
#relative error variance
PGY3_c <- read.csv("PGY3_c.csv")
rel_err_var <- PGY3_c[19,4]/n
#absolute error variance
abs_err_var <- PGY3_c[19,3]/n +PGY3_c[19,4]/n
#calculate generalization coefficient
gen_coef <- PGY3_c[19,2]/(PGY3_c[19,2]+ rel_err_var)
#calculate dependability coefficient
dep_coef <- PGY3_c[19,2]/(PGY3_c[19,2]+ abs_err_var)
#calculate SEM
SEM <- sqrt(abs_err_var)

label <- "Delta"
data <- data.frame(n=n,SEM=SEM) %>%
  mutate(label=label)


p_PGY3 <- ggplot(data, aes(n=n,y=SEM)) +
  geom_text(aes(n,SEM,label=label), parse=TRUE) +
  scale_y_continuous(breaks = c(0,0.2,0.4,0.6,0.8,1,1.2,1.4,1.6))+
  xlab("Number of Procedures") + 
  ylab("Standard Error of Measurement") + 
  labs(title = "PGY 3", tag = "b") +
  theme_bw()+
  theme(axis.line = element_line(colour = "black"),
        panel.background = element_blank(),
        panel.border = element_blank(),
        legend.title = element_blank())+
  theme(plot.title = element_text(hjust = 0.5,size=10),
        plot.tag = element_text())

############################ PGY4 #########################################
n <- 1:20
#relative error variance
PGY4_c <- read.csv("PGY4_c.csv")
rel_err_var <- PGY4_c[19,4]/n
#absolute error variance
abs_err_var <- PGY4_c[19,3]/n +PGY4_c[19,4]/n
#calculate generalization coefficient
gen_coef <- PGY4_c[19,2]/(PGY4_c[19,2]+ rel_err_var)
#calculate dependability coefficient
dep_coef <- PGY4_c[19,2]/(PGY4_c[19,2]+ abs_err_var)
#calculate SEM
SEM <- sqrt(abs_err_var)

label <- "Delta"
data <- data.frame(n=n,SEM=SEM) %>%
  mutate(label=label)


p_PGY4 <- ggplot(data, aes(n=n,y=SEM)) +
  geom_text(aes(n,SEM,label=label), parse=TRUE) +
  scale_y_continuous(breaks = c(0,0.2,0.4,0.6,0.8,1,1.2,1.4,1.6))+
  xlab("Number of Procedures") + 
  ylab("Standard Error of Measurement") + 
  labs(title = "PGY 4", tag = "c") +
  theme_bw()+
  theme(axis.line = element_line(colour = "black"),
        panel.background = element_blank(),
        panel.border = element_blank(),
        legend.title = element_blank())+
  theme(plot.title = element_text(hjust = 0.5,size=10),
        plot.tag = element_text())

############################ PGY5 #########################################
n <- 1:20
#relative error variance
PGY5_c <- read.csv("PGY5_c.csv")
rel_err_var <- PGY5_c[33,4]/n
#absolute error variance
abs_err_var <- PGY5_c[33,3]/n +PGY5_c[33,4]/n
#calculate generalization coefficient
gen_coef <- PGY5_c[33,2]/(PGY5_c[33,2]+ rel_err_var)
#calculate dependability coefficient
dep_coef <- PGY5_c[33,2]/(PGY5_c[33,2]+ abs_err_var)
#calculate SEM
SEM <- sqrt(abs_err_var)

label <- "Delta"

data <- data.frame(n=n,SEM=SEM) %>%
  mutate(label=label)


p_PGY5 <- ggplot(data, aes(n=n,y=SEM)) +
  geom_text(aes(n,SEM,label=label), parse=TRUE) +
  scale_y_continuous(breaks = c(0,0.2,0.4,0.6,0.8,1,1.2,1.4,1.6))+
  xlab("Number of Procedures") + 
  ylab("Standard Error of Measurement") + 
  labs(title = "PGY 5", tag = "d") +
  theme_bw()+
  theme(axis.line = element_line(colour = "black"),
        panel.background = element_blank(),
        panel.border = element_blank(),
        legend.title = element_blank())+
  theme(plot.title = element_text(hjust = 0.5,size=10),
        plot.tag = element_text())


tiff("SEM.tiff", unit="in",width=7, height=5, res=500)
ggarrange(p_PGY2,p_PGY3,p_PGY4,p_PGY5)
dev.off()